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Abstract 



The Perdew-Burke-Ernzerhof generalized gradient approximation to the 
density functional theory is tested with respect to sensitivity to the choice of 
the value of the parameter k, which is associated to the degree of localization of 
the exchange-correlation hole. A study of structural and dynamical properties 
of four selected ferroelectric perovskites is presented. The originally proposed 
value of ac=0.804 works well for some solids, whereas for the ABO3 perovskites 
it must be decreased in order to predict equilibrium lattice parameters in 
good agreement with experiments. The effects on the structural instabilities 
and zone center phonon modes are examined. The need of varying k from 
one system to another reflects the fact that the localization of the exchange- 
correlation hole is system dependent, and the sensitivity of the structural 
properties to its actual value illustrates the necessity of finding a universal 
function for k. 



INTRODUCTION 



Due to their wide variety of physical properties, the ABO3 perovskites constitute a 
particularly interesting family for studies of structural phase transitions. These depend 
strongly on the composition of the compounds, as for instance on the chemical nature of 
the A and B ions in pure samples, as well as on the relative composition in the case of 
solid solutions. Furthermore, the temperature driven structural phase transitions of a given 
pure system are strongly affected by the crystal volume and thus by application of external 
pressure. 

The electronic, structural and dynamical properties of ferroelectric perovskites have been 
the subject of intense theoretical work, but only recently has first principles density func- 
tional theory been applied with success to analyze them at a microscopic level. The local 
(spin-) density approximation, L(S)DA, has provided a useful framework to elucidate im- 
portant aspects of the underlying physics of these oxides, providing important quantitative 
information about electronic charge distributions, character of electronic bondings, crys- 
tal structure, phonons, structural instabilities, etc. For example, frozen phonon and linear 
response approaches have been applied to the study of the lattice dynamics of many per- 
ovskites, such as KNb0 3 |H8, BaTi0 3 HI, KTa0 3 00], SrTi0 3 §, etc ; providing 
considerable insight into the nature of the soft-mode total energy surface. 

The main drawback of the LDA when applied to perovskites is the systematic under- 
estimation of the equilibrium volume, and several properties are poorly reproduced at this 
volume. In fact, in prototypical materials, like KNb0 3 and BaTi0 3 , ferroelectricity is 
strongly suppressed at the LDA equilibrium volume. Therefore, L(S)DA studies often use 
a simple ad hoc correction consisting in using the experimental rather than the calculated 
equilibrium volumes. This is just one example of shortcomings that have motivated a num- 
ber of attempts to go beyond the LSDA in the search of a better framework. The application 
of the Perdew and Wang (PW91) version || of the Generalized Gradient Approximation 
(GGA) functional to KNb0 3 and BaTi0 3 showed that the equilibrium volumes of the two 
compounds are overcorrected; although in KNb0 3 GGA yielded an equilibrium volume 



that is much closer to the experimental volume than the LSDA result, in BaTiOs the GGA 
error was almost as large as found within the LDA. |TIJ More recently, the Weighted Density 
Approximation (WDA) has been applied to KNbO"3 yielding an equilibrium volume in close 
agreement with experiment. 

In the Kohn-Sham density functional theory only the exchange-correlation energy 
Exc =Ex+Ec which is a functional of the electron spin densities must be approximated, 
and for slowly varying densities, n, it can be expressed as the volume integrals of n times 
e™ % * in the LSDA case and /( n f , n j ,V n | ,V n J, ) for the GGA case. For practical 
calculations the exchange-correlation energy density of a uniform electron gas, e""^ ( n f , 
n I), and / must be parametrized. The form of e"™^ is now well established but which is 
the best choice of / is still under debate. 

In the recent work of Perdew, Burke and Ernzerhof (PBE) ]12[ a simplified form of the 
GGA-/ was presented which improves over the previous PW91 in an accurate description of 
the linear response of the uniform gas, behaving correctly under uniform scaling and giving a 
smoother potential. The new PBE retains the correct features of LSDA and combines them 
with the most energetically important features of gradient corrected nonlocality. Studies of 



atomization energies for small molecules |L2[ gave essentially the same results as PW91. 

However, a remainig uncertainty relates to the value of the parameter k, in the en- 
hancement factor Fx(s) which is directly associated to the degree of localization of the 
exchange- correlation hole. (Here the variable s is a measure of the relative density gradient, 
s = \Wn\/2kpn kp giving the Fermi wavenumber of an electron gas of density n). In their 
original work PBE proposed Fx(s) = 1 + k — k/(1 + fis 2 /^) which satisfied the inequality 
F x < 1.804 with k = 0.804 and with the value of [i = 0.21951. Zhang and Yang found 
13| that the results for several atoms and molecules were improved by increasing k beyond 



the originally proposed value, 0.804. But, this does not hold for all type of bonds |L4 



and it may well happen that applications to solids appear to be improved when smaller 
values of k are used. The parameter k might be a weak function of the reduced Laplacian, 
K=g(V 2 n/(2kp) 2 n), which would also satisfy the condition 'D' in the PBE paper of uniform 
density scaling, i.e. that Ex scales like A. 



A study of LSDA and PBE-GGA (k=0.804) has recently been made for s — p materials 
n|, 3d and Ad transition metals [T5| and another, similar to the present, also including 



magnetic effects, has been made for the 3d, 4d and 5d transition metal series jnj . 

In the present work we take four different compounds as representative examples of 
the ferroelectric perovskites in their cubic phases. The PBE proposal of the GGA for the 
exchange correlation energy within DFT is tested against its intrinsic uncertainty in choosing 
the coefficient k and also compared to the LSDA results. For this we have selected three 
structure related properties which relate to the physics of the ferroelectric instabilities and 
phase transitions: (1) the energy-volume curves from which lattice parameters and bulk 
moduli are derived, (2) r 15 phonon modes for the cubic structure, and (3) the tendency 
to suffer a ferroelectric transition when the atoms are displaced according to the soft-mode 
displacement pattern. 

The selected perovskites are KNb03 , BaTi03 , SrTi03 and KTaOs , since they exem- 
plify the rich variety of structure related and dynamical effects of the perovkite materials. 
These oxides have been extensively studied because of their possible technological applica- 
tions, some of which are related to their ferroelectric properties. But a basic understanding 
of the mechanisms of their various transformations of phase and structure is also of great 
interest . 

The two compounds, KNbOs and BaTiOs , exhibit the same sequence of phase transi- 
tions upon cooling from the high-temperature cubic paraelectric phase to slightly distorted 
ferroelectric structures with tetragonal, orthorhombic, and rhombohedral symmetry. Their 
transition temperatures have an approximately constant ratio : T(KNb03)/T(BaTiOs) = 
1.6. Thus the dynamics of these two materials should be similar but with different energy 
scales. SrTiOs has the simple cubic structure at high temperatures going through an an- 
tiferrodistortive transition at 104 K to a tetragonal phase in which the oxygen octahedra 
are rotated in opposite senses in neighboring unit cells. KTaOs remains cubic perovskite 
down to low temperatures. These two last materials present a softening of a transverse optic 
phonon which correspond to the ferroelectric mode, but its frequency remains finite at all 
temperatures. This fact induced these two compounds to be termed incipient ferroelectrics. 



COMPUTATIONAL DETAILS 



The calculations presented in this work were performed within the LDA and the PBE- 
GGA to density functional theory, using the full-potential Linear Augmented Plane Wave 
(LAPW) method. In this method no shape approximation on either the potential or the 



electronic charge density is made. We use the WIEN97 implementation of the method fl8 



which allows the inclusion of local orbitals (LO) in the basis, improving upon linearization 
and making possible a consistent treatment of semicore and valence states in one energy 
window hence ensuring proper orthogonality ||19|| . Ceperley and Alder [2(J L(S)DA exchange- 



correlation contribution, as parametrized by Perdew and Zunger |21|, was used. 

The muffin-tin sphere radii (Ri) = 2.0, 2.0, 2.0 ,1.95, 1.90, 1.90 and 1.50 a.u. were used for 
K, Ba, Sr, Ti, Nb, Ta, and O, respectively. The value of the parameter RK max , which controls 
the size of the basis sets in these calculations, was chosen to be 8 for all systems studied. 
This gives well converged basis sets consisting of approximately 1200 LAPW functions plus 
local orbitals for the different systems. In order to obtain sufficient accuracy with inclusion 
of the GGA functionals a relatively high plane wave cut-off energy (256 Ry) was necessary 
in order to obtain a converged interstitial representation of the potential. This choice of 
parameters was justified by performing calculations for other Ri values and by increasing 
RKmax- Although an increase of RK max from 8 to 9 had no significant effects on most of 
physical properties that we studied, convergence of the cohesive energy to a precision of 
1 mRy per formula unit required RK max — 10 with the choice of Ri given above. 

We introduced LO to include the following orbitals in the basis set: K-3s and 3p, Ba-5s, 
bp and Ad, Sr-As and Ap, Ti-3s and 3p, Nb-4s and Ap, Ta-5s, bp and Af and 0-2s. 

Integrations in reciprocal space were performed using the special points method. We used 
6x6x6 meshes which represent 250 /c-points in the first Brillouin zone. This corresponds 
to 28 special k-points in the irreducible wedge for the rhombohedral structure, 18 for the 
tetragonal and 10 in the cubic. Convergence tests indicate that only small changes result 
from increasing to a denser k-mesh. 

Of the two allowed symmetries for T point phonon modes in the cubic perovskite structure 



we only studied the triply degenerate Tis modes. There are three of these, not counting the 
zero frequency acoustic mode. These are infrared active and include the "ferroelectric soft" 
mode. We determined the phonon frequencies and polarizations for this particular symmetry 
by calculating atomic forces for several small displacements ( ~ 0.01 A ) consistent with 
the symmetry and small enough to be in the linear regime. From the force as a function of 
displacement the dynamical matrix was constructed and diagonalized. 



RESULTS 

Table I lists the lattice parameters and bulk moduli as derived from energy- volume curves 
obtained within the LDA as well as the PBE-GGA (with k = 0.804 as originally proposed by 



PBE [0), together with the experimental values. The equilibrium volumes obtained from 
the LDA calculations are 4.3%, 5.0 % , 2.3 % and 2.0 % smaller than experiment for KNb0 3 
, BaTi0 3 , SrTi0 3 and KTa0 3 , respectively. This agrees with the generally observed 
tendency to overbinding in this approximation. The LDA bulk moduli are overestimated 
when evaluated at these too small volumes. The values calculated at the experimental 
volumes (see the values in brackets) agree better with the experiments. 

As is the general trend in the application of GGA functionals to solids, it expands bonds, 
an effect that sometimes corrects and sometimes overcorrects the LSDA prediction. All GGA 
calculations with k=0.804 overestimate the equilibrium volumes, namely by 1.6 % , 1.6 % , 
4.2 % and 4.0 % for KNb0 3 , BaTi0 3 , SrTi0 3 and KTa0 3 , respectively. We agree with 



previous calculations by Singh [TI| for KNb0 3 , but our results for BaTi0 3 are different. 
Indeed, we obtained for both oxides that the GGA yields to an equilibrium volume that 
is much closer to the experimental volume than the LSDA result. However, for the two 
incipient ferroelectrics the GGA errors are twice as large as the LSDA ones. We have not 
found any GGA calculations for SrTi0 3 and KTa0 3 to compare with. 

In Fig. 1 the values of V/V G (V G is the corresponding experimental value) are given 
for the four perovskites when the value of k in the PBE-GGA functional is decreased from 
0.804 to 0.3. As can be seen both BaTi0 3 and KNb0 3 would give perfect determination 



of the lattice parameter (V/V G =1) for k 0.6. In the case of SrTi0 3 and KTa0 3 the 
value should be further reduced to m 0.4. The need of varying k from one system to another 
reflects the fact that the localization of the exchange-correlation hole is system dependent. 
Nevertheless, it is important to investigate the sensitivity of the vibrational properties and 
energetics which are involved in the ferroelectric instabilities by the application of the PBE- 
GGA functional and the reduction of the k value. This constitutes a hard test on the quality 
of any ab-initio scheme, since the curvature of the total energy surface over the manifold of 
varios atomic displacements is much more sensitive to the details of the calculations than 
merely the position of the total energy minimun. 

To clarify the effect of the GGA functional on the phonon energies, we have performed 
frozen phonon calculations for the four selected perovskites at their experimental lattice con- 
stants, and examined the effects of choosing different exchange-correlation approximations, 
i.e. LDA and PBE-GGA with different k values. The calculated frequencies of the r 15 modes 
and the experimental values are shown in Table II. Two conclusions can be drawn from these 
calculations. Firstly, the GGA functional hardens the phonon frequencies (as compared to 
the LDA frequencies) of the four perovskites. This hardening produces a slight reduction of 
the errors, since the LSDA provides phonon frequencies which are understimated by 10-20% 
compared with experiments. The second conclusion concerns the parameter k. It is evident 
from Table II that the effect introduced on the GGA phonon frequencies by the modification 
of k is negligible. 

This picture also holds for the low-frequency T 15 mode, i.e. the ferroelectric mode. In 
the cases of KNb0 3 and KTa0 3 , the GGA produces a small change of its frequency when 
it is compared with the LSDA result. However, a remarkable hardening is observed in the 
titanates where the GGA increases the imaginary frequency in BaTi0 3 by ps 40 cm -1 and 
stabilises the imaginary frequency mode in SrTi0 3 . In the latter material, the potential 
well where atoms move is very flat so the small differences in energetics between LSDA and 
PBE-GGA might give different predictions which would modify the underlying physics. 

As a by-product of the frozen phonon calculation, we show in Table III the calculated 
eigenvectors of KNb0 3 to exemplify the effect produced on the phonon polarizations by the 



application of the GGA. Note that the eigenvectors are practically unchanged. The same 
occurs for the other three perovskites. 

Finally, to test the sensitivity of the energetics involved in the ferroelectric instabilities 
when the different exchange-correlation functionals are used, we performed total energy 
calculations as a function of the ferroelectric soft mode displacement pattern. In Figure 2, 
we show the energy as a function of the amplitude of (111) ferroelectric mode displacements 
for KNb0 3 and BaTi0 3 . The calculations were performed at the experimental lattice 
constants. In the case of KNb0 3 , the LSDA and GGA (with /t=0.804 and 0.6) yield a 
clear ferroelectric instability with similar energetics and displacements, and with an energy 
lowering of 1.8 mRy/cell. However, for BaTi0 3 the LSDA yields an energy lowering 
of 1.75 mRy/cell while the GGA yields 1.25 mRy/cell, which is consistent with the above 
mentioned hardening of the soft-mode frequency produced by the GGA. The energetics in 
both materials are very little dependent on the k value used in the GGA functional. Although 
the transition temperatures are not associated directly with the energy barriers shown in 
Figure 2, it is interesting to note that while for the LSDA the energy barrier ratio of the 
two compounds is E(KNb0 3 )/E(BaTi0 3 ) « 1, for the GGA E(KNb0 3 )/E(BaTi0 3 ) « 1.5 
which is much closer to the transition temperatures ratio T(KNb0 3 )/T(BaTi0 3 ) = 1.6. 

CONCLUSIONS 

Results of LSDA and PBE-GGA calculations have been reported for structural and 
dynamical properties relevant for the physics of the ferroelectric instabilities of perovskites. 
KNb0 3 , BaTi0 3 , SrTi0 3 and KTa0 3 were taken as representative example materials. In 
the case of PBE-GGA calculations were made with different values of the coefficient k which 
is related to the localization of the exchange-correlation hole. The usual underestimation of 
LSDA for the equilibrium volume is obtained (an average of approximately 3% for the four 
materials) in complete accord with all other theoretical studies. The PBE-GGA, on the other 
hand, overestimates the equilibrium volumes when /t=0.804 (the originally proposed value) 
is used. For k values of 0.6 in the cases of KNb0 3 and BaTi0 3 and 0.4 in SrTi0 3 and 



KTa0 3 a correct equilibrium volume (V/V o =1.0) can be obtained. The need of varying k 
from one system to another reflects the fact that the localisation of the exchange-correlation 
hole is system dependent. However, the vibrational properties and energetics involved in 
the ferroelectric instabilities do not depend on the use of different k values. This fact would 
permit an ad hoc correction of the theoretical equilibrium volume by selecting the adequate 
value of k for each particular system. However, this procedure would not lead to a fully 
ab-initio, free-paremeters, calculation scheme, indicating the necessity of finding a universal 
function for /t=g(V 2 n/(2kp) 2 n). 
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FIGURES 



Figure 1: Equilibrium volumes relative to the experimental ones of KNb0 3 , BaTi0 3 , 
SrTiO"3 and KTa0 3 as a function of the PBE-GGA intrinsic parameter k. 

Figure 2: Total energy as a function of ferroelectric displacement along the (111) direction 
for KNb0 3 and BaTi0 3 , calculated at the experimental lattice constant for different 
exchange-correlation energies: LSDA and PBE-GGA (k = 0.804 and 0.6). 



TABLES 

TABLE I. Lattice parameter (in a.u.) and Bulk modulus (in GPa) obtained for LSDA and 
PBE-GGA(k = 0.804) and compared with experimental values. In parantheses: bulk modulus at 
the experimental equilibrium volume. 

LSDA PBE-GGA (k = 0.804) Experiment 



KNbO. 



a 
B 



7.48 
206(155) 



7.63 
171(186) 



7.59 a 
138 b 



BaTiO. 



a 
B 



7.44 
195(155) 



7.61 
160(173) 



7.57 a 



SrTiO, 



a 
B 



7.30 
204(176) 



7.46 
167(194) 



7.358 a 
179 c 



KTaO ? 



a 
B 



7.475 
222(192) 



7.625 
183(213) 



7.526 a 
220 d 



a corresponds to the experimental equilibrium volume extrapolated to T=0. 
b calculated from cubic elastic constants, Ref. |22 . 



c Ref. [|3|. 
d Ref. m. 



TABLE II. Frequencies of the modes obtained for LSDA and PBE-GGA with two differents 
values of k. n e quii. is the value of k for which the calculated equilibrium volume agrees with 
experiment (see Fig. 1.) 



Compound 




Frequency (cm 1 ) 






LDA 


GGA (k=0.804) GGA (K=K eguiL ) 


Expt. 



KNb0 3 2 Hi 197i 195i soft 

166 182 179 198 a 

466 478 478 521 a 



BaTi0 3 239i 203i 205i soft 

163 168 167 182 b 

454 463 461 482 b 



SrTi0 3 57i 64 52 soft 

157 166 166 175 c 

532 551 545 545 c 



KTa0 3 99 107 109 soft 

172 189 181 197 d 

529 550 545 546 d 



a Infrared reflectivity measurements at 710 K, Ref ||25|| . 



Infrared reflectivity measurements at 420 K, Ref [26]. 
c Infrared reflectivity measurements at room temperature, Ref. [27]. 
d Hyper Raman scattering at room temperature, Ref. 



TABLE III. Frequencies and eigenvectors of the T 15 modes in KNDO3 using LSDA and 



PBE-GGA (k=0.804 and 0.6). 





Freq. (cm 1 ) 




Ei£ 


;envectors 








K 


Nb 


01 


02 and 03 


KNb0 3 


210i 


0.015 


-0.58 


0.6 


0.39 


LSDA, a=7.59au. 


166 


-0.88 


0.36 


0.146 


0.178 




466 


0.027 


-0.11 


-0.73 


0.48 


KNbOg 


196i 


0.042 


-0.59 


0.57 


0.40 


GGA, a=7.59au. 


182 


-0.88 


0.35 


0.178 


0.18 


k = 0.804 


478 


0.008 


-0.08 


-0.75 


0.47 


KNb0 3 


195i 


0.035 


-0.59 


0.567 


0.40 


GGA, a=7.59au. 


179 


-0.88 


0.35 


0.17 


0.18 


k = 0.6 


478 


0.012 


-0.083 


-0.75 


0.46 



1.04 - 



1.03 - 



1.02 - 



1.01 - 



1.00 - 



0.99 




0.3 0.4 0.5 0.6 0.7 0.8 



Kappa 



— BaTiO, 



LDA 

GGA K=0.8 
GGA K=0.6 




0.00 0.05 0.00 0.05 0.10 

B displacement relative to A (A) 



